// hnumcolamd.h
#pragma once

#ifndef __HNUMCOLAMD_H__
#define __HNUMCOLAMD_H__

#include "hnumdef.h"
#include "hwinexception.h"

namespace harlinn
{
    /*
     * This is a C++ implementation of the "column approximate minimum degree ordering" algorthm.
     *
     * This implementation is based on the original c implementation by
     * Stefan I. Larimore and Timothy A. Davis
     *
     *      colamd:  An approximate minimum degree column ordering algorithm.
     *
     *      Purpose:
     *      
	 *      Colamd computes a permutation Q such that the Cholesky factorization of
	 *      (AQ)'(AQ) has less fill-in and requires fewer floating point operations
	 *      than A'A.  This also provides a good ordering for sparse partial
	 *      pivoting methods, P(AQ) = LU, where Q is computed prior to numerical
	 *      factorization, and P is computed during numerical factorization via
	 *      conventional partial pivoting with row interchanges.  Colamd is the
	 *      column ordering method used in SuperLU, part of the ScaLAPACK library.
	 *      It is also available as user-contributed software for Matlab 5.2,
	 *      available from MathWorks, Inc. (http://www.mathworks.com).  This
	 *      routine can be used in place of COLMMD in Matlab.  By default, the \
	 *      and / operators in Matlab perform a column ordering (using COLMMD)
	 *      prior to LU factorization using sparse partial pivoting, in the
	 *      built-in Matlab LU(A) routine.
     *      
     *      Authors:
     *      
	 *      The authors of the code itself are Stefan I. Larimore and Timothy A.
	 *      Davis (davis@cise.ufl.edu), University of Florida.  The algorithm was
	 *      developed in collaboration with John Gilbert, Xerox PARC, and Esmond
	 *      Ng, Oak Ridge National Laboratory.
     *      
     *      Date:
     *      
	 *      August 3, 1998.  Version 1.0.
     *      
     *      Acknowledgements:
     *      
	 *      This work was supported by the National Science Foundation, under
	 *      grants DMS-9504974 and DMS-9803599.
     *      
     *      Notice:
     *      
	 *      Copyright (c) 1998 by the University of Florida.  All Rights Reserved.
     *      
	 *      THIS MATERIAL IS PROVIDED AS IS, WITH ABSOLUTELY NO WARRANTY
	 *      EXPRESSED OR IMPLIED.  ANY USE IS AT YOUR OWN RISK.
     *      
	 *      Permission is hereby granted to use or copy this program for any
	 *      purpose, provided the above notices are retained on all copies.
	 *      User documentation of any code that uses this code must cite the
	 *      Authors, the Copyright, and "Used by permission."  If this code is
	 *      accessible from within Matlab, then typing "help colamd" or "colamd"
	 *      (with no arguments) must cite the Authors.  Permission to modify the
	 *      code and to distribute modified code is granted, provided the above
	 *      notices are retained, and a notice that the code was modified is
	 *      included with the above copyright notice.  You must also retain the
	 *      Availability information below, of the original version.
     *      
	 *      This software is provided free of charge.
     *      
     *      Availability:
     *      
	 *      This file is located at
     *      
	 *      	http://www.cise.ufl.edu/~davis/colamd/colamd.c
     *      
	 *      The colamd.h file is required, located in the same directory.
	 *      The colamdmex.c file provides a Matlab interface for colamd.
	 *      The symamdmex.c file provides a Matlab interface for symamd, which is
	 *      a symmetric ordering based on this code, colamd.c.  All codes are
	 *      purely ANSI C compliant (they use no Unix-specific routines, include
	 *      files, etc.).
     *
     */
    namespace numerics
    {
        using namespace harlinn::windows;

        class ColumnApproximateMinimumDegreeOrdering
        {
            

            static const int STATS = 20;
            static const int DENSE_ROW = 0;
            static const int DENSE_COL = 1;
            static const int DEFRAG_COUNT = 2;
            static const int JUMBLED_COLS = 3;

            static const int EMPTY	= -1;

            // rowData & column status
            static const int ALIVE	= 0;
            static const int DEAD	= UINT32_MAX;

            // Column status
            static const int DEAD_PRINCIPAL		= UINT32_MAX;
            static const int DEAD_NON_PRINCIPAL	= UINT32_MAX - 1;
            static const int FLAG = 0x80000000;
            static const int MASK = 0x7FFFFFFF;

            
            


            inline static void Debug(const char*,...)
            {

            }


            struct ColumnData
            {
                // index for workArea of first row in this column, or DEAD 
			    // if column is dead 
                int start;
                // number of rows in this column 
                int length;
                union
                {
                    // number of original columns represented by this column, if the column is alive 
	                int thickness;
                    // parent in parent tree super-column structure, if the column is dead 
	                int parent;
                };
                union
                {
                    // the score used to maintain heap, if col is alive 
	                int score;	
                    // pivot ordering of this column, if col is dead 
	                int order;	
                };
                union
                {
                    // head of a hash bucket, if col is at the head of a degree list 
	                int headhash;
                    // hash value, if col is not in a degree list 
	                int hash;
                    // previous column in degree list, if col is in a degree list (but not at the head of a degree list) 
	                int prev;	
                };
                union
                {
                    // next column, if col is in a degree list
	                int degree_next;
                    // next column, if col is in a hash list
	                int hash_next;
                };

                bool IsDead() const
                {
                    return start >= ColumnApproximateMinimumDegreeOrdering::DEAD_NON_PRINCIPAL;
                }

                bool IsDeadPrincipal() const
                {
                    return start == ColumnApproximateMinimumDegreeOrdering::DEAD_PRINCIPAL;
                }

                bool IsAlive() const
                {
                    return start < ColumnApproximateMinimumDegreeOrdering::DEAD_NON_PRINCIPAL;
                }

                void KillPrincipal() { start = ColumnApproximateMinimumDegreeOrdering::DEAD_PRINCIPAL; }
                void KillNonPrincipal()	{ start = ColumnApproximateMinimumDegreeOrdering::DEAD_NON_PRINCIPAL; }

            };

            struct RowData
            {
                // index for workArea of first col in this row 
                int start;		
                // number of principal columns in this row 
                int length;
                union
                {
                    // number of principal & non-principal columns in row 
	                int degree;	
                    // used as a row pointer in init_rows_cols () 
	                int p;
                };
                union
                {
                    // for computing set differences and marking dead rows
	                int mark;
                    // first column in row (used in garbage collection)
	                int first_column;
                };

                bool IsDead() const
                {
                    return mark == ColumnApproximateMinimumDegreeOrdering::DEAD;
                }

                bool IsAlive() const
                {
                    return mark < ColumnApproximateMinimumDegreeOrdering::DEAD;
                }

                void Kill() { mark = ColumnApproximateMinimumDegreeOrdering::DEAD; }

            };


            int numberOfRows;
            int numberOfColumns;
            int workAreaLength;
            int* workArea;
            ColumnData* columnData;
            RowData* rowData;
            double denseRow; // 0.5
            double denseCol; // 0.5

        public:
            ColumnApproximateMinimumDegreeOrdering()
                : numberOfRows(0),
                  numberOfColumns(0),
                  workAreaLength(0),
                  workArea(nullptr),
                  columnData(nullptr),
                  rowData(nullptr),
                  denseRow(0.5),
                  denseCol(0.5)
            {}
        private:


            // InitializeRowsAndColumns
            // 
            // Takes the column form of the matrix in workArea and creates the row form of the
            // matrix.  Also, row and column attributes are stored in the columnData and rowData
            // structs.  If the columns are un-sorted or contain duplicate row indexes,
            // this routine will also sort and remove duplicate row indexes from the
            // column form of the matrix.  
            //
            // int columnIndexes []	array of the column indexes for the workArea, of size numberOfColumns+1
            // 
            // returns true if the columns are jumbled 

            bool InitializeRowsAndColumns(int columnIndexes[])
            {
                int *cp;			    // a column pointer
                int *cp_end;		    // a pointer to the end of a column
                int *rp;			    // a row pointer
                int *rp_end;		    // a pointer to the end of a row
                int last_start;	// start index of previous column in workArea
                int start;			// start index of column in workArea
                bool jumbled_columns;   // indicates if columns are jumbled

                // Initialize the columnData 
                last_start = 0 ;
                for (int i = 0 ; i < numberOfColumns ; i++)
                {
                    start = columnIndexes[i];
                    // verify that the column indexes are ordered
                    if (start < last_start)
                    {
                        throw ArgumentException(Exception::Format("Column pointers must not be decreasing, last %d current %d",last_start,start));
                    }
                    columnData[i].start = start;
                    columnData[i].length = columnIndexes[i+1] - start;
                    columnData[i].thickness = 1;
                    columnData[i].score = 0;
                    columnData[i].prev = EMPTY;
                    columnData[i].degree_next = EMPTY;
                    last_start = start;
                }
                // verify the column index for the last column
                if (columnIndexes[numberOfColumns] < last_start)
                {
                    throw ArgumentException(Exception::Format("Column pointers must not be decreasing, last %d current %d",columnIndexes[numberOfColumns],last_start));
                }

                

                jumbled_columns = false;

                for (int i = 0 ; i < numberOfRows ; i++)
                {
                    rowData[i].length = 0;
                    rowData[i].mark = -1;
                }

                // Scan columns, compute row degrees, and check row indexes

                for (int i = 0 ; i < numberOfColumns ; i++)
                {
                    int previousRowIndex = -1 ;

                    cp = &workArea [columnIndexes [i]];
                    cp_end = &workArea [columnIndexes [i+1]];

                    while (cp < cp_end)
                    {
                        int rowIndex = *cp;

                        // make sure row indexes within range
                        if (rowIndex < 0 || rowIndex >= numberOfRows)
                        {
                            throw ArgumentException(Exception::Format("row index out of range, column %d row %d previous row %d",i, rowIndex, previousRowIndex));
                        }
                        else if (rowIndex <= previousRowIndex )
                        {
                            // row indexes are not sorted or repeated, thus cols are jumbled
                            jumbled_columns = true;
                        }
                        // prevent repeated row from being counted 
                        if (rowData [rowIndex].mark != i)
                        {
                            rowData [rowIndex].length++ ;
                            rowData [rowIndex].mark = i;
                            previousRowIndex = rowIndex;
                        }
                        else
                        {
                            // this is a repeated entry in the column, it will be removed 
                            columnData[i].length--;
                        }
                        cp++;
                    }
                }

                // Compute row pointers

                // row form of the matrix starts directly after the column form of matrix in workArea 
                rowData [0].start = columnIndexes [numberOfColumns];
                rowData [0].p = rowData [0].start ;
                rowData [0].mark = -1 ;
                for (int i = 1 ; i < numberOfRows ; i++)
                {
                    rowData [i].start = rowData [i-1].start + rowData [i-1].length;
                    rowData [i].p = rowData [i].start;
                    rowData [i].mark = -1;
                }

                // Create row form

                if (jumbled_columns)
                {
                    // if cols jumbled, watch for repeated row indexes
                    for (int i = 0 ; i < numberOfColumns ; i++)
                    {
                        cp = &workArea [columnIndexes[i]];
                        cp_end = &workArea [columnIndexes[i+1]];
                        while (cp < cp_end)
                        {
                            int row = *cp++ ;
                            if (rowData [row].mark != i)
                            {
                                workArea [(rowData [row].p)++] = i;
                                rowData [row].mark = i;
                            }
                        }
                    }
                }
                else
                {
                    // if cols not jumbled, we don't need the mark (this is faster) 
                    for (int i = 0 ; i < numberOfColumns ; i++)
                    {
                        cp = &workArea [columnIndexes[i]];
                        cp_end = &workArea [columnIndexes[i+1]];
                        while (cp < cp_end)
                        {
                            workArea [(rowData [*cp++].p)++] = i;
                        }
                    }
                }

                // Clear the row marks and set row degrees

                for (int i = 0 ; i < numberOfRows ; i++)
                {
                    rowData [i].mark = 0 ;
                    rowData [i].degree = rowData [i].length;
                }

                // See if we need to re-create columns

                if (jumbled_columns)
                {

                    // Compute col pointers 

                    // col form of the matrix starts at workArea [0].
                    // Note, we may have a gap between the col form and the row
                    // form if there were duplicate entries, if so, it will be
                    // removed upon the first garbage collection 
                    columnData [0].start = 0 ;
                    columnIndexes[0] = columnData [0].start ;
                    for (int i = 1 ; i < numberOfColumns ; i++)
                    {
                        // note that the lengths here are for pruned columns, i.e. 
                        // no duplicate row indexes will exist for these columns 
                        columnData [i].start = columnData [i-1].start + columnData [i-1].length ;
                        columnIndexes[i] = columnData [i].start ;
                    }

                    // Re-create col form

                    for (int i = 0 ; i < numberOfRows; i++)
                    {
                        rp = &workArea [rowData [i].start];
                        rp_end = rp + rowData [i].length ;
                        while (rp < rp_end)
                        {
                            workArea [(columnIndexes[*rp++])++] = i;
                        }
                    }
                    return true;
                }
                else
                {
                    // no columns jumbled (this is faster)
                    return false;
                }
            }


            // InitializeScoring
            // 
            //
            // Kills dense or empty columns and rows, calculates an initial score for
            // each column, and places all columns in the degree lists.
            // 
            // Parameters:
            // 
            //  int head [],		of size numberOfColumns+1
            //  int *p_n_row2,		number of non-dense, non-empty rows
            //  int *p_n_col2,		number of non-dense, non-empty columns
            //  int *p_max_deg		maximum row degree

            void InitializeScoring(int head [],int *p_n_row2,int *p_n_col2,int *p_max_deg)
            {
                int c;			// a column index 
                int r, row;	// a row index 
                int *cp;		    // a column pointer 
                int *cp_end;	    // a pointer to the end of a column
                int *new_cp;	    // new column pointer
                int col_length;// length of pruned column
                int score;		    // current column score
                int n_col2;	// number of non-dense, non-empty columns
                int n_row2;	// number of non-dense, non-empty rows
                int dense_row_count;// remove rows with more entries than this
                int dense_col_count;// remove cols with more entries than this
                int min_score;		// smallest column score
                int max_deg;		// maximum row degree
                int next_col;		// Used to add to degree list.
#ifndef NDEBUG
                int debug_count ;	// debug only.
#endif

                dense_row_count = max (0, min (int(denseRow * numberOfColumns), numberOfColumns));
                dense_col_count = max (0, min (int(denseCol * numberOfRows), numberOfRows));
                max_deg = 0;
                n_col2 = numberOfColumns;
                n_row2 = numberOfRows;

                // Kill empty columns

                // Put the empty columns at the end in their natural, so that LU 
                // factorization can proceed as far as possible. 
                for (c = numberOfColumns-1 ; c >= 0 ; c--)
                {
                    if (columnData [c].length == 0)
                    {
                        // this is a empty column, kill and order it last
                        n_col2--;
                        columnData[c].order = n_col2;
                        columnData[c].KillPrincipal();
                    }
                }

                // Kill dense columns

                // Put the dense columns at the end, in their natural order
                for (c = numberOfColumns-1 ; c >= 0 ; c--)
                {
                    // skip any dead columns
                    if (columnData[c].IsDead())
                    {
                        continue ;
                    }
                    int deg = columnData [c].length ;
                    if (deg > dense_col_count)
                    {
                        // this is a dense column, kill and order it last 
                        columnData [c].order = --n_col2 ;
                        // decrement the row degrees 
                        cp = &workArea [columnData [c].start] ;
                        cp_end = cp + columnData [c].length ;
                        while (cp < cp_end)
                        {
                            rowData [*cp++].degree--;
                        }
                        columnData[c].KillPrincipal();
                    }
                }
                Debug("Dense and null columns killed: %d\n", numberOfColumns - n_col2);

                // Kill dense and empty rows

                for (r = 0 ; r < numberOfRows ; r++)
                {
                    int deg = rowData [r].degree ;
                    assert (deg >= 0 && deg <= numberOfColumns) ;
                    if (deg > dense_row_count || deg == 0)
                    {
                        // kill a dense or empty row 
                        rowData [r].Kill();
                        --n_row2 ;
                    }
                    else
                    {
                        // keep track of max degree of remaining rows
                        max_deg = max (max_deg, deg) ;
                    }
                }
                Debug("Dense and null rows killed: %d\n", numberOfRows - n_row2);

                // Compute initial column scores

                // At this point the row degrees are accurate.  They reflect the number
                // of "live" (non-dense) columns in each row.  No empty rows exist.
                // Some "live" columns may contain only dead rows, however.  These are
                // pruned in the code below.

                // now find the initial matlab score for each column
                for (c = numberOfColumns-1 ; c >= 0 ; c--)
                {
                    // skip dead column
                    if (columnData[c].IsDead())
                    {
                        continue ;
                    }
                    score = 0;
                    cp = &workArea [columnData [c].start];
                    new_cp = cp;
                    cp_end = cp + columnData [c].length;
                    while (cp < cp_end)
                    {
                        // get a row
                        row = *cp++ ;
                        // skip if dead
                        if (rowData[row].IsDead())
                        {
                            continue;
                        }
                        // compact the column
                        *new_cp++ = row ;
                        // add row's external degree
                        score += rowData [row].degree - 1 ;
                        // guard against integer overflow 
                        score = min (score, numberOfColumns) ;
                    }
                    // determine pruned column length
                    col_length = (int) (new_cp - &workArea [columnData [c].start]) ;
                    if (col_length == 0)
                    {
                        // a newly-made null column (all rows in this col are "dense" 
                        // and have already been killed)
                        Debug("Newly null killed: %d\n", c);
                        columnData [c].order = --n_col2 ;
                        columnData[c].KillPrincipal();
                        
                    }
                    else
                    {
                        // set column length and set score
                        assert (score >= 0) ;
                        assert (score <= numberOfColumns) ;
                        columnData [c].length = col_length ;
                        columnData [c].score = score ;
                    }
                }
                Debug("Dense, null, and newly-null columns killed: %d\n",numberOfColumns-n_col2);

                // At this point, all empty rows and columns are dead.  All live columns 
                // are "clean" (containing no dead rows) and simplicial (no supercolumns
                // yet).  Rows may contain dead columns, but all live rows contain at
                // least one live column.

                // Initialize degree lists

#ifndef NDEBUG
                debug_count = 0 ;
#endif

                // clear the hash buckets
                for (c = 0 ; c <= numberOfColumns ; c++)
                {
                    head [c] = EMPTY;
                }
                min_score = numberOfColumns ;
                // place in reverse order, so low column indexes are at the front
                // of the lists.  This will encourage natural tie-breaking
                for (c = numberOfColumns-1 ; c >= 0 ; c--)
                {
                    // only add principal columns to degree lists
                    if (columnData[c].IsAlive())
                    {
                        Debug("place %d score %d minscore %d ncol %d\n", c, columnData [c].score, min_score, numberOfColumns);

                        // Add columns score to DList

                        score = columnData [c].score ;

                        assert (min_score >= 0) ;
                        assert (min_score <= numberOfColumns) ;
                        assert (score >= 0) ;
                        assert (score <= numberOfColumns) ;
                        assert (head [score] >= EMPTY) ;

                        // now add this column to dList at proper score location 
                        next_col = head [score] ;
                        columnData [c].prev = EMPTY ;
                        columnData [c].degree_next = next_col ;

                        // if there already was a column with the same score, set its
                        // previous pointer to this new column 
                        if (next_col != EMPTY)
                        {
                            columnData [next_col].prev = c ;
                        }
                        head [score] = c ;

                        // see if this score is less than current min 
                        min_score = min (min_score, score) ;

#ifndef NDEBUG
                        debug_count++ ;
#endif
                    }
                }


                // Return number of remaining columns, and max row degree 

                *p_n_col2 = n_col2 ;
                *p_n_row2 = n_row2 ;
                *p_max_deg = max_deg ;
            }


            // DetermineOrder
            //
            // Order the principal columns of the supercolumn form of the matrix
            // (no supercolumns on input).  Uses a minimum approximate column minimum
            // degree ordering method.  Not user-callable.
            //
            // Parameters
            // 
            //  int head [],		of size numberOfColumns+1 
            //  int n_col2,			Remaining columns to order 
            //  int max_deg,		Maximum row degree 
            //  int pfree			index of first free slot (2*nnz on entry) 
            //  
            // returns the number of garbage collections 

            int DetermineOrder(int head [],int n_col2,int max_deg,int pfree)
            {
                int k;			// current pivot ordering step 
                int pivot_col;		// current pivot column 
                int *cp;			// a column pointer 
                int *rp;			// a row pointer 
                int pivot_row;		// current pivot row 
                int *new_cp;		// modified column pointer 
                int *new_rp;		// modified row pointer 
                int pivot_row_start;	// pointer to start of pivot row 
                int pivot_row_degree;	// # of columns in pivot row 
                int pivot_row_length;	// # of supercolumns in pivot row 
                int pivot_col_score;	// score of pivot column 
                int needed_memory;		// free space needed for pivot row 
                int *cp_end;		// pointer to the end of a column 
                int *rp_end;		// pointer to the end of a row 
                int row;			// a row index 
                int col;			// a column index 
                int max_score;		// maximum possible score 
                int cur_score;		// score of current column 
                int hash;		    // hash value for supernode detection 
                int head_column;	// head of hash bucket 
                int first_col;		// first column in hash bucket 
                int tag_mark;		// marker value for mark array 
                int row_mark;		// rowData [row].mark 
                int set_difference;	// set difference size of row with pivot row 
                int min_score;		// smallest column score 
                int col_thickness;		// "thickness" (# of columns in a supercol) 
                int max_mark;		// maximum value of tag_mark
                int pivot_col_thickness;	// number of columns represented by pivot col
                int prev_col;		// Used by Dlist operations.
                int next_col;		// Used by Dlist operations.
                int ngarbage;		// number of garbage collections performed
#ifndef NDEBUG
                int debug_d ;		// debug loop counter
                int debug_step = 0 ;	// debug loop counter
#endif

                // Initialization and clear mark
                max_mark = INT_MAX - numberOfColumns ;	// INT_MAX defined in <limits.h> 
                tag_mark = ClearMark();
                min_score = 0 ;
                ngarbage = 0 ;
                Debug("Ordering.. n_col2=%d\n", n_col2);

                // Order the columns

                for (k = 0 ; k < n_col2 ; )// 'k' is incremented below 
                {
                    // Select pivot column, and order it

                    // make sure degree list isn't empty 
                    assert (min_score >= 0) ;
                    assert (min_score <= numberOfColumns) ;
                    assert (head [min_score] >= EMPTY) ;

#ifndef NDEBUG
                    for (debug_d = 0 ; debug_d < min_score ; debug_d++)
                    {
                        assert (head [debug_d] == EMPTY) ;
                    }
#endif

                    // get pivot column from head of minimum degree list 
                    while (head [min_score] == EMPTY && min_score < numberOfColumns)
                    {
                        min_score++ ;
                    }
                    pivot_col = head [min_score] ;
                    assert (pivot_col >= 0 && pivot_col <= numberOfColumns) ;
                    next_col = columnData [pivot_col].degree_next ;
                    head [min_score] = next_col ;
                    if (next_col != EMPTY)
                    {
                        columnData [next_col].prev = EMPTY ;
                    }

                    assert (columnData[pivot_col].IsAlive()) ;
                    Debug("Pivot col: %d\n", pivot_col);

                    // remember score for defrag check 
                    pivot_col_score = columnData [pivot_col].score ;

                    // the pivot column is the kth column in the pivot order 
                    columnData [pivot_col].order = k ;

                    // increment order count by column thickness 
                    pivot_col_thickness = columnData [pivot_col].thickness ;
                    k += pivot_col_thickness ;
                    assert (pivot_col_thickness > 0) ;

                    // Garbage_collection, if necessary

                    needed_memory = min (pivot_col_score, numberOfColumns - k) ;
                    if (pfree + needed_memory >= workAreaLength)
                    {
                        pfree = CompactWorkArea(&workArea [pfree]);
                        ngarbage++ ;
                        // after garbage collection we will have enough 
                        assert (pfree + needed_memory < workAreaLength) ;
                        // garbage collection has wiped out the rowData[].mark array
                        tag_mark = ClearMark() ;
                    }

                    // Compute pivot row pattern

                    // get starting location for this new merged row 
                    pivot_row_start = pfree ;

                    // initialize new row counts to zero 
                    pivot_row_degree = 0 ;

                    // tag pivot column as having been visited so it isn't included 
                    // in merged pivot row 
                    columnData [pivot_col].thickness = -pivot_col_thickness ;

                    // pivot row is the union of all rows in the pivot column pattern 
                    cp = &workArea [columnData [pivot_col].start] ;
                    cp_end = cp + columnData [pivot_col].length ;
                    while (cp < cp_end)
                    {
                        // get a row 
                        row = *cp++ ;
                        Debug("Pivot col pattern %d %d\n", rowData[row].IsAlive(), row);
                        // skip if row is dead 
                        if (rowData[row].IsDead())
                        {
                            continue ;
                        }
                        rp = &workArea [rowData [row].start] ;
                        rp_end = rp + rowData [row].length ;
                        while (rp < rp_end)
                        {
                            // get a column 
                            col = *rp++ ;
                            // add the column, if alive and untagged 
                            col_thickness = columnData [col].thickness ;
                            if (col_thickness > 0 && columnData[col].IsAlive() )
                            {
                                // tag column in pivot row 
                                columnData [col].thickness = -col_thickness ;
                                assert (pfree < workAreaLength) ;
                                // place column in pivot row 
                                workArea [pfree++] = col ;
                                pivot_row_degree += col_thickness ;
                            }
                        }
                    }

                    // clear tag on pivot column 
                    columnData [pivot_col].thickness = pivot_col_thickness ;
                    max_deg = max (max_deg, pivot_row_degree) ;

                    // Kill all rows used to construct pivot row 

                    // also kill pivot row, temporarily 
                    cp = &workArea [columnData [pivot_col].start] ;
                    cp_end = cp + columnData [pivot_col].length ;
                    while (cp < cp_end)
                    {
                        // may be killing an already dead row 
                        row = *cp++ ;
                        Debug("Kill row in pivot col: %d\n", row);
                        rowData [row].Kill();
                    }

                    // Select a row index to use as the new pivot row 

                    pivot_row_length = pfree - pivot_row_start ;
                    if (pivot_row_length > 0)
                    {
                        // pick the "pivot" row arbitrarily (first row in col) 
                        pivot_row = workArea [columnData [pivot_col].start] ;
                        Debug("Pivotal row is %d\n", pivot_row);
                    }
                    else
                    {
                        // there is no pivot row, since it is of zero length 
                        pivot_row = EMPTY ;
                        assert (pivot_row_length == 0) ;
                    }
                    assert (columnData [pivot_col].length > 0 || pivot_row_length == 0) ;

                    // Approximate degree computation 

                    // Here begins the computation of the approximate degree.  The column 
                    // score is the sum of the pivot row "length", plus the size of the 
                    // set differences of each row in the column minus the pattern of the 
                    // pivot row itself.  The column ("thickness") itself is also 
                    // excluded from the column score (we thus use an approximate 
                    // external degree).

                    // The time taken by the following code (compute set differences, and
                    // add them up) is proportional to the size of the data structure
                    // being scanned - that is, the sum of the sizes of each column in
                    // the pivot row.  Thus, the amortized time to compute a column score
                    // is proportional to the size of that column (where size, in this
                    // context, is the column "length", or the number of row indexes
                    // in that column).  The number of row indexes in a column is
                    // monotonically non-decreasing, from the length of the original
                    // column on input to colamd.

                    // Compute set differences

                    Debug("** Computing set differences phase. **\n");

                    // pivot row is currently dead - it will be revived later. 

                    Debug("Pivot row: ");
                    // for each column in pivot row 
                    rp = &workArea [pivot_row_start] ;
                    rp_end = rp + pivot_row_length ;
                    while (rp < rp_end)
                    {
                        col = *rp++ ;
                        assert (columnData[col].IsAlive() && col != pivot_col) ;
                        Debug("columnData: %d\n", col);

                        // clear tags used to construct pivot row pattern 
                        col_thickness = -columnData [col].thickness ;
                        assert (col_thickness > 0) ;
                        columnData [col].thickness = col_thickness ;

                        // Remove column from degree list

                        cur_score = columnData [col].score ;
                        prev_col = columnData [col].prev ;
                        next_col = columnData [col].degree_next ;
                        assert (cur_score >= 0) ;
                        assert (cur_score <= numberOfColumns) ;
                        assert (cur_score >= EMPTY) ;
                        if (prev_col == EMPTY)
                        {
                            head [cur_score] = next_col ;
                        }
                        else
                        {
                            columnData [prev_col].degree_next = next_col ;
                        }
                        if (next_col != EMPTY)
                        {
                            columnData [next_col].prev = prev_col ;
                        }

                        // Scan the column 

                        cp = &workArea [columnData [col].start] ;
                        cp_end = cp + columnData [col].length ;
                        while (cp < cp_end)
                        {
                            /* get a row */
                            row = *cp++ ;
                            
                            /* skip if dead */
                            if (rowData [row].IsDead())
                            {
                                continue ;
                            }
                            row_mark = rowData [row].mark ;
                            assert (row != pivot_row) ;
                            set_difference = row_mark - tag_mark ;
                            // check if the row has been seen yet
                            if (set_difference < 0)
                            {
                                assert (rowData [row].degree <= max_deg) ;
                                set_difference = rowData [row].degree ;
                            }
                            // subtract column thickness from this row's set difference 
                            set_difference -= col_thickness ;
                            assert (set_difference >= 0) ;
                            // absorb this row if the set difference becomes zero 
                            if (set_difference == 0)
                            {
                                Debug("aggressive absorption. rowData: %d\n", row);
                                rowData [row].Kill();
                            }
                            else
                            {
                                // save the new mark
                                rowData [row].mark = set_difference + tag_mark ;
                            }
                        }
                    }


                    // Add up set differences for each column

                    Debug("** Adding set differences phase. **\n");

                    // for each column in pivot row 
                    rp = &workArea [pivot_row_start] ;
                    rp_end = rp + pivot_row_length ;
                    while (rp < rp_end)
                    {
                        // get a column 
                        col = *rp++ ;
                        assert (columnData[col].IsAlive() && col != pivot_col) ;
                        hash = 0 ;
                        cur_score = 0 ;
                        cp = &workArea [columnData [col].start] ;
                        // compact the column 
                        new_cp = cp ;
                        cp_end = cp + columnData [col].length ;

                        Debug("Adding set diffs for columnData: %d.\n", col);

                        while (cp < cp_end)
                        {
                            // get a row 
                            row = *cp++ ;
                            assert(row >= 0 && row < numberOfRows) ;
                            
                            // skip if dead 
                            if (rowData [row].IsDead())
                            {
                                continue ;
                            }
                            row_mark = rowData [row].mark ;
                            assert (row_mark > tag_mark) ;
                            // compact the column 
                            *new_cp++ = row ;
                            // compute hash function 
                            hash += row ;
                            // add set difference 
                            cur_score += row_mark - tag_mark ;
                            // integer overflow... 
                            cur_score = min (cur_score, numberOfColumns) ;
                        }

                        // recompute the column's length 
                        columnData [col].length = (int) (new_cp - &workArea [columnData [col].start]) ;

                        // Further mass elimination

                        if (columnData [col].length == 0)
                        {
                            Debug("further mass elimination. columnData: %d\n", col);
                            // nothing left but the pivot row in this column 
                            columnData[col].KillPrincipal();
                            pivot_row_degree -= columnData [col].thickness;
                            assert (pivot_row_degree >= 0) ;
                            // order it 
                            columnData [col].order = k ;
                            // increment order count by column thickness 
                            k += columnData [col].thickness ;
                        }
                        else
                        {
                            // Prepare for supercolumn detection

                            Debug("Preparing supercol detection for columnData: %d.\n", col);

                            // save score so far
                            columnData [col].score = cur_score ;

                            // add column to hash table, for supercolumn detection 
                            hash %= numberOfColumns + 1 ;

                            Debug(" Hash = %d, numberOfColumns = %d.\n", hash, numberOfColumns);
                            assert (hash <= numberOfColumns) ;

                            head_column = head [hash] ;
                            if (head_column > EMPTY)
                            {
                                // degree list "hash" is non-empty, use prev of 
                                // first column in degree list as head of hash bucket 
                                first_col = columnData [head_column].headhash ;
                                columnData [head_column].headhash = col ;
                            }
                            else
                            {
                                // degree list "hash" is empty, use head as hash bucket 
                                first_col = - (head_column + 2) ;
                                head [hash] = - (col + 2) ;
                            }
                            columnData [col].hash_next = first_col ;

                            // save hash function in columnData [col].hash 
                            columnData [col].hash = (int) hash ;
                            assert (columnData[col].IsAlive()) ;
                        }
                    }

                    // The approximate external column degree is now computed. 

                    // Supercolumn detection 

                    Debug("** Supercolumn detection phase. **\n");

                    detect_super_cols (head, pivot_row_start, pivot_row_length);

                    // Kill the pivotal column

                    columnData[pivot_col].KillPrincipal();

                    // Clear mark

                    tag_mark += (max_deg + 1) ;
                    if (tag_mark >= max_mark)
                    {
                        Debug("clearing tag_mark\n");
                        tag_mark = ClearMark() ;
                    }

                    // Finalize the new pivot row, and column scores 

                    Debug("** Finalize scores phase. **\n");

                    // for each column in pivot row
                    rp = &workArea [pivot_row_start] ;
                    // compact the pivot row 
                    new_rp = rp ;
                    rp_end = rp + pivot_row_length ;
                    while (rp < rp_end)
                    {
                        col = *rp++ ;
                        // skip dead columns
                        if (columnData[col].IsDead())
                        {
                            continue;
                        }
                        *new_rp++ = col ;
                        // add new pivot row to column 
                        workArea [columnData [col].start + (columnData [col].length++)] = pivot_row;

                        // retrieve score so far and add on pivot row's degree.
                        // (we wait until here for this in case the pivot 
                        // row's degree was reduced due to mass elimination).
                        cur_score = columnData [col].score + pivot_row_degree;

                        // calculate the max possible score as the number of
                        // external columns minus the 'k' value minus the
                        // columns thickness
                        max_score = numberOfColumns - k - columnData [col].thickness;

                        // make the score the external degree of the union-of-rows 
                        cur_score -= columnData [col].thickness;

                        // make sure score is less or equal than the max score
                        cur_score = min (cur_score, max_score);
                        assert (cur_score >= 0);

                        // store updated score 
                        columnData [col].score = cur_score;

                        // Place column back in degree list

                        assert (min_score >= 0);
                        assert (min_score <= numberOfColumns);
                        assert (cur_score >= 0);
                        assert (cur_score <= numberOfColumns);
                        assert (head [cur_score] >= EMPTY);
                        next_col = head [cur_score];
                        columnData [col].degree_next = next_col;
                        columnData [col].prev = EMPTY;
                        if (next_col != EMPTY)
                        {
                            columnData [next_col].prev = col;
                        }
                        head [cur_score] = col;

                        // see if this score is less than current min 
                        min_score = min (min_score, cur_score) ;

                    }

                    // Resurrect the new pivot row

                    if (pivot_row_degree > 0)
                    {
                        // update pivot row length to reflect any cols that were killed 
                        // during super-col detection and mass elimination 
                        rowData [pivot_row].start  = pivot_row_start ;
                        rowData [pivot_row].length = (int) (new_rp - &workArea[pivot_row_start]) ;
                        rowData [pivot_row].degree = pivot_row_degree ;
                        rowData [pivot_row].mark = 0 ;
                        // pivot row is no longer dead 
                    }
                }

                // All principal columns have now been ordered

                return ngarbage;
            }


            // order_children
            // 
            // The find_ordering routine has ordered all of the principal columns (the
            // representatives of the supercolumns).  The non-principal columns have not
            // yet been ordered.  This routine orders those columns by walking up the
            // parent tree (a column is a child of the column which absorbed it).  The
            // final permutation vector is then placed in p [0 ... numberOfColumns-1], with p [0]
            // being the first column, and p [numberOfColumns-1] being the last.  It doesn't look
            // like it at first glance, but be assured that this routine takes time linear
            // in the number of columns.  Although not immediately obvious, the time
            // taken by this routine is O (numberOfColumns), that is, linear in the number of
            // columns.
            //
            // Parameters
            // 
            //  int numberOfColumns,			number of columns of workArea 
            //  ColumnData columnData [],		of size numberOfColumns+1 
            //  int p []			p [0 ... numberOfColumns-1] is the column permutation

            void order_children(int p [])
            {
                int i;			// loop counter for all columns
                int c;			// column index 
                int parent;		// index of column's parent 
                int order;		// column's order 

                // Order each non-principal column

                for (i = 0 ; i < numberOfColumns ; i++)
                {
                    // find an un-ordered non-principal column 
                    assert (columnData[i].IsDead()) ;
                    if (!columnData[i].IsDeadPrincipal() && columnData [i].order == EMPTY)
                    {
                        parent = i ;
                        // once found, find its principal parent 
                        do
                        {
                            parent = columnData [parent].parent ;
                        } while (!columnData[parent].IsDeadPrincipal()) ;

                        // now, order all un-ordered non-principal columns along path 
                        // to this parent.  collapse tree at the same time 
                        c = i ;
                        // get order of parent 
                        order = columnData [parent].order ;

                        do
                        {
                            assert (columnData [c].order == EMPTY) ;

                            // order this column 
                            columnData [c].order = order++ ;
                            // collaps tree 
                            columnData [c].parent = parent ;

                            // get immediate parent of this column 
                            c = columnData [c].parent ;

                            // continue until we hit an ordered column.  There are 
                            // guarranteed not to be anymore unordered columns 
                            // above an ordered column 
                        } while (columnData [c].order == EMPTY) ;

                        // re-order the super_col parent to largest order for this group 
                        columnData [parent].order = order ;
                    }
                }

                // Generate the permutation
                for (c = 0 ; c < numberOfColumns ; c++)
                {
                    p [columnData [c].order] = c ;
                }
            }


            // detect_super_cols
            // 
            //  Detects supercolumns by finding matches between columns in the hash buckets.
            //  Check amongst columns in the set workArea [row_start ... row_start + row_length-1].
            //  The columns under consideration are currently *not* in the degree lists,
            //  and have already been placed in the hash buckets.
            //  
            //  The hash bucket for columns whose hash function is equal to h is stored
            //  as follows:
            //  
            //  if head [h] is >= 0, then head [h] contains a degree list, so:
            //  
            //  head [h] is the first column in degree bucket h.
            //  columnData [head [h]].headhash gives the first column in hash bucket h.
            //  
            //  otherwise, the degree list is empty, and:
            //  
            //  -(head [h] + 2) is the first column in hash bucket h.
            //  
            //  For a column c in a hash bucket, columnData [c].prev is NOT a "previous
            //  column" pointer.  columnData [c].hash is used instead as the hash number
            //  for that column.  The value of columnData [c].hash_next is the next column
            //  in the same hash bucket.
            //  
            //  Assuming no, or "few" hash collisions, the time taken by this routine is
            //  linear in the sum of the sizes (lengths) of each column whose score has
            //  just been computed in the approximate degree computation.
            //  
            // Parameters
            // 
            //  int head []		    head of degree lists and hash buckets
            //  int row_start		pointer to set of columns to check
            //  int row_length		number of columns to check
            //

            void detect_super_cols(int head [],int row_start,int row_length)
            {
                int hash;		// hash # for a column 
                int *rp;		// pointer to a row 
                int c;			// a column index 
                int super_c;	// column index of the column to absorb into 
                int *cp1;		// column pointer for column super_c
                int *cp2;		// column pointer for column c
                int length;		// length of column super_c
                int prev_c;		// column preceding c in hash bucket
                int i;			// loop counter
                int *rp_end;	// pointer to the end of the row
                int col;		// a column index in the row to check
                int head_column;// first column in hash bucket or degree list
                int first_col;	// first column in hash bucket

                // Consider each column in the row

                rp = &workArea [row_start];
                rp_end = rp + row_length;
                while (rp < rp_end)
                {
                    col = *rp++;
                    if (columnData[col].IsDead())
                    {
                        continue;
                    }

                    // get hash number for this column
                    hash = columnData[col].hash;
                    assert (hash <= numberOfColumns);

                    // Get the first column in this hash bucket

                    head_column = head [hash];
                    if (head_column > EMPTY)
                    {
                        first_col = columnData [head_column].headhash;
                    }
                    else
                    {
                        first_col = - (head_column + 2);
                    }

                    // Consider each column in the hash bucket

                    for (super_c = first_col ; super_c != EMPTY; super_c = columnData [super_c].hash_next)
                    {
                        assert (columnData[super_c].IsAlive() ) ;
                        assert (columnData [super_c].hash == hash) ;
                        length = columnData [super_c].length ;

                        // prev_c is the column preceding column c in the hash bucket 
                        prev_c = super_c ;

                        // Compare super_c with all columns after it 

                        for (c = columnData [super_c].hash_next ;c != EMPTY ; c = columnData [c].hash_next)
                        {
                            assert (c != super_c);
                            assert (columnData[c].IsAlive());
                            assert (columnData [c].hash == hash);

                            // not identical if lengths or scores are different
                            if (columnData [c].length != length ||columnData [c].score != columnData [super_c].score)
                            {
                                prev_c = c;
                                continue;
                            }

                            // compare the two columns
                            cp1 = &workArea [columnData [super_c].start];
                            cp2 = &workArea [columnData [c].start];

                            for (i = 0 ; i < length ; i++)
                            {
                                // the columns are "clean" (no dead rows)
                                assert (rowData[(*cp1)].IsAlive() )  ;
                                assert (rowData[(*cp2)].IsAlive() )  ;
                                // row indexes will same order for both supercols, 
                                // no gather scatter nessasary 
                                if (*cp1++ != *cp2++)
                                {
                                    break ;
                                }
                            }

                            // the two columns are different if the for-loop "broke" 
                            if (i != length)
                            {
                                prev_c = c ;
                                continue ;
                            }

                            // Got it!  two columns are identical 

                            assert (columnData [c].score == columnData [super_c].score) ;

                            columnData [super_c].thickness += columnData [c].thickness ;
                            columnData [c].parent = super_c ;
                            columnData[c].KillNonPrincipal();
                            // order c later, in order_children()
                            columnData [c].order = EMPTY ;
                            // remove c from hash bucket 
                            columnData [prev_c].hash_next = columnData [c].hash_next ;
                        }
                    }

                    // Empty this hash bucket

                    if (head_column > EMPTY)
                    {
                        // corresponding degree list "hash" is not empty 
                        columnData [head_column].headhash = EMPTY ;
                    }
                    else
                    {
                        // corresponding degree list "hash" is empty 
                        head [hash] = EMPTY ;
                    }
                }
            }


            // CompactWorkArea
            //
            // Defragments and compacts columns and rows in the workspace workArea.  Used when
            // all avaliable memory has been used while performing row merging.  Returns
            // the index of the first free position in workArea, after garbage collection.  The
            // time taken by this routine is linear is the size of the array workArea, which is
            // itself linear in the number of nonzeros in the input matrix.
            // 
            // Parameters
            // 
            //     int *pfree			&workArea [0] ... pfree is in use
            //
            // returns the new value of pfree 
            //
            int CompactWorkArea(int *pfree)
            {
                int *psrc;		 // source pointer 
                int *pdest;		 // destination pointer 
                int j;		 // counter 
                int r;		 // a row index 
                int c;		 // a column index 
                int length; // length of a row or column 

#ifndef NDEBUG
                int debug_rows ;
                Debug("Defrag..\n");
                for (psrc = &workArea[0] ; psrc < pfree ; psrc++) assert (*psrc >= 0) ;
                debug_rows = 0 ;
#endif
                // Defragment the columns
                pdest = &workArea[0];
                for (c = 0 ; c < numberOfColumns ; c++)
                {
                    if (columnData[c].IsAlive())
                    {
                        psrc = &workArea[columnData[c].start];

                        // move and compact the column 
                        assert (pdest <= psrc) ;
                        columnData[c].start = (int) (pdest - &workArea[0]);
                        length = columnData[c].length;
                        for (j = 0 ; j < length ; j++)
                        {
                            r = *psrc++ ;
                            if (rowData[r].IsAlive() )
                            {
                                *pdest++ = r ;
                            }
                        }
                        columnData[c].length = (int)(pdest - &workArea [columnData[c].start]);
                    }
                }

                // Prepare to defragment the rows
                for (r = 0 ; r < numberOfRows ; r++)
                {
                    if (rowData [r].IsAlive())
                    {
                        if (rowData [r].length == 0)
                        {
                            // this row is of zero length.  cannot compact it, so kill it 
                            Debug("Defrag row kill\n");
                            rowData [r].Kill();
                        }
                        else
                        {
                            // save first column index in rowData [r].first_column 
                            psrc = &workArea [rowData [r].start];
                            rowData [r].first_column = *psrc;
                            assert (rowData[r].IsAlive());
                            // flag the start of the row
                            *psrc |= FLAG;
#ifndef NDEBUG
                            debug_rows++;
#endif
                        }
                    }
                }

                // Defragment the rows

                psrc = pdest ;
                while (psrc < pfree)
                {
                    // find a negative number ... the start of a row 
                    if (*psrc & FLAG)
                    {
                        // get back the row index by removing the flag bit
                        r = *psrc & MASK;
                        assert (r >= 0 && r < numberOfRows) ;
                        // restore first column index 
                        *psrc = rowData [r].first_column ;
                        assert (rowData[r].IsAlive()) ;

                        // move and compact the row 
                        assert (pdest <= psrc) ;
                        rowData [r].start = (int) (pdest - &workArea [0]) ;
                        length = rowData [r].length ;
                        for (j = 0 ; j < length ; j++)
                        {
                            c = *psrc++ ;
                            if ( columnData[c].IsAlive() )
                            {
                                *pdest++ = c ;
                            }
                        }
                        rowData [r].length = (int) (pdest - &workArea [rowData [r].start]) ;
#ifndef NDEBUG
                        debug_rows-- ;
#endif
                    }
                    psrc++;
                }
                // ensure we found all the rows 
                assert (debug_rows == 0) ;

                // Return the new value of pfree

                return ((int) (pdest - &workArea [0])) ;
            }


            // ClearMark 
            //
            //  Clears the rowData [].mark array, and returns the new tag_mark.
            //  Return value is the new tag_mark
            //  
            // Parameters
            // 
            //     int numberOfRows		number of rows in workArea
            //     RowData rowData []	rowData [0 ... numberOfRows-1].mark is set to zero 
            //
            // return the new value for tag_mark 

            int ClearMark()
            {
                int r;

                Debug("Clear mark\n");
                for (r = 0 ; r < numberOfRows ; r++)
                {
                    if (rowData[r].IsAlive())
                    {
                        rowData [r].mark = 0 ;
                    }
                }
                return 1;
            }


        public:
            static size_t Recommended(size_t nnz, size_t numberOfRows, size_t numberOfColumns)
            {
                size_t minimum; // bare minimum requirements 
                size_t recommended;	// recommended value of Alen 

                minimum =
	                2 * (nnz)		// for workArea 
	                + (((numberOfColumns) + 1) * sizeof (ColumnData) / sizeof (int))	// for columnData 
	                + (((numberOfRows) + 1) * sizeof (RowData) / sizeof (int))	// for rowData 
	                + numberOfColumns			// minimum elbow room to guarrantee success 
	                + STATS ;	// for output statistics 

                /* recommended is equal to the minumum plus enough memory to keep the */
                /* number garbage collections low */
                recommended = minimum + nnz/5;

                return recommended;
            }


            // Parameters
            // 
            //  int theNumberOfRows		number of rows in theMatrixData
            //  int theNumberOfColumns	    number of columns in theMatrixData
            //  int theMatrixDataLength 	length of theMatrixData
            //  int theMatrixData[] 			row indexes of theMatrixData
            //  int p []		            column indexes in theMatrixData
            //
            // returns true if successful 

            bool Execute(int theNumberOfRows,int theNumberOfColumns,int theMatrixDataLength,int theMatrixData [],int p [])
            {
                int nnz;			// nonzeros in theMatrixData
                int Row_size;		// size of rowData [], in integers
                int Col_size;		// size of columnData [], in integers
                int elbow_room;	// remaining free space
                int n_col2;		// number of non-dense, non-empty columns
                int n_row2;		// number of non-dense, non-empty rows
                int ngarbage;		// number of garbage collections performed
                int max_deg;		    // maximum row degree
                int init_result;	    // return code from initialization

#ifndef NDEBUG
                debug_colamd = 0; // no debug printing 
                // get "D" environment variable, which gives the debug printing level 
                if (getenv ("D")) debug_colamd = atoi (getenv ("D")) ;
                Debug("debug version, D = %d (THIS WILL BE SLOOOOW!)\n", debug_colamd);
#endif

                // Check the input arguments

                if (!theNumberOfRows || !theNumberOfColumns || !theMatrixDataLength || !theMatrixData || !p)
                {
                    // numberOfRows and numberOfColumns must be non-negative, workArea and p must be present 
                    Debug("error! %d %d %d\n", theNumberOfRows, theNumberOfColumns, theMatrixDataLength);
                    return (FALSE) ;
                }

                this->workArea = theMatrixData;
                numberOfRows = theNumberOfRows;
                numberOfColumns = theNumberOfColumns;
                workAreaLength = theMatrixDataLength;

                nnz = p [numberOfColumns] ;
                if (nnz < 0 || p [0] != 0)
                {
                    // nnz must be non-negative, and p [0] must be zero 
                    Debug("colamd error! %d %d\n", nnz, p [0]);
                    return false;
                }


                // Allocate the rowData and columnData arrays from array workArea
                Col_size = (numberOfColumns + 1) * sizeof (ColumnData) / sizeof (int);
                Row_size = (numberOfRows + 1) * sizeof (RowData) / sizeof (int);
                elbow_room = theMatrixDataLength - (2*nnz + Col_size + Row_size);

                if (elbow_room < numberOfColumns + STATS)
                {
                    // not enough space in array workArea to perform the ordering
                    Debug("error! elbow_room %d, %d\n", elbow_room,numberOfColumns);
                    return (FALSE) ;
                }
                workAreaLength = 2*nnz + elbow_room;
                columnData  = (ColumnData *) &workArea[workAreaLength];
                rowData  = (RowData *) &workArea[workAreaLength + Col_size];

                // Construct the row and column data structures

                init_result = InitializeRowsAndColumns(p);
                if (init_result == -1)
                {
                    // input matrix is invalid
                    Debug("error! matrix invalid\n");
                    return (FALSE) ;
                }

                // Initialize scores, kill dense rows/columns

                InitializeScoring(p, &n_row2, &n_col2, &max_deg) ;

                // Order the supercolumns

                ngarbage = DetermineOrder(p,n_col2, max_deg, 2*nnz);

                // Order the non-principal columns

                order_children (p);

                // Return statistics in theMatrixData

                for (int i = 0 ; i < STATS ; i++)
                {
                    theMatrixData [i] = 0 ;
                }
                theMatrixData [DENSE_ROW] = numberOfRows - n_row2 ;
                theMatrixData [DENSE_COL] = numberOfColumns - n_col2 ;
                theMatrixData [DEFRAG_COUNT] = ngarbage ;
                theMatrixData [JUMBLED_COLS] = init_result ;

                return true;
            }

        };


    };

};



#endif //__HNUMCOLAMD_H__
